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We simulated the field-dependent magnetization m{H,T) and the uniform susceptibility x{H,T) 
of classical Heisenberg antiferromagnets in the chain and square-lattice geometry using Monte 
Carlo methods. The results confirm the singular behavior of x{H,T) at small T,H: 
YimT^o\imH~.ox{H,T) = l/(2Jo)(l - 1/D) and lim^^o limr^o x(^^, = l/(2Jo), where D = ?, 
is the number of spin components, Jo = zJ , and z is the number of nearest neighbors. A good 
agreement is achieved in a wide range of temperatures T and magnetic fields H with the first-order 
1/D expansion resuhs [D. A. Garanin, J. Stat. Phys. 83, 907 (1996)]. 
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In recent years, investigations of two-dimensional an- 
tiferromagnets concentrated primarily on the quantum 
model with S — 1/2. A practical reason for that is its 
possible relevance for the high-temperature superconduc- 
tivity. On the other hand, the identification with the 
quantum nonlinear sigma model (QNLtrM) in the low- 
energy sector allowed using field-theory methods 
Although the QNLctM resuhs for the 5" = 1/2 model 
proved to be in a good argeement with quantum Monte 
Carlo (QMC) simulations (see, e.g., Ref. ||]), the re- 
quirement of low energies confines the validity region 
of the QNLcrM to rather low temperatures already for 
S' > 1. High-temperature series expansions (HTSE) for 
5 > 1 @ and QMC simulations @ for = 1 in the 
experimentally relevant temperature range, as well as 
experiments on model substances with 1 < 5 < 5/2, 
showed much better accord with the pure-quantum self- 
consistent harmonic approximation (PQSCHA) than 
with the field-theoretical QNLcrM predictions. In con- 
trast to the QNLcrM, the PQSCHA maps a quantum sys- 
tem on the corresponding classical system on the lattice, 
which, in turn, can be studied by classical MC simula- 
tions or other methods. The parameters of these classical 
Hamiltonians are renormalized by quantum fluctuations 
and given by explicit analytical expressions. 

The above arguments show that in most cases the clas- 
sical model can be used as a good starting point for study- 
ing quantum systems. In fact, most of nontrivial features 
of two-dimensional antiferromagnets, such as impossibil- 
ity of ordering at nonzero temperatures in the isotropic 
case, are universal and appear already at the classical 
level. The main theoretical problem is that due to Gold- 
stone modes, a simple spin-wave theory at T ^ JS^ is 
inapplicable to two-dimensional magnets. 

Despite their importance, classical antiferromagnets 
received much less attention than the quantum S* = 1/2 



model. In particular, the initial uniform susceptibility 
x(r) for the square lattice having a flat maximum at 
T ^ J has been simulated for S* = 1/2 in Refs. |^,^ and 
for S* = 1 in Ref. but there are no results for the 
classical model yet! For the latter, only the old MC data 
for the energy |^ are available up to now. 

On the other hand, classical magnets can be theoreti- 
cally studied with the help of the 1/D expansion, where 
D is the number of spin components . In Ref. , 

x{T) has been calculated for the square lattice and lin- 
ear chain to first order in 1/D for all temperatures, the 
solution interpolating between the exact result at T = 
and the leading terms of the HTSE at high tempera- 
tures. In contrast, the low-energy approaches such as 
"Schwinger-boson mean-field theory" p3[ | or "modified 
spin-wave theory" ||l^ break down at T > J and fail 
to reproduce the maximum of x{T)- It should be noted 
that for quantum magnets there is a method consisting in 
the expansion in powers of 1/7V where N is the number 
of fiavors in the Schwinger-boson technique ||l^. This 
method, which is nonequivalent to the 1/D expansion in 
the limit 5' — *■ oo, is supposed to work for all T, in con- 
trast to the low-energy QNLctM. Unfortunately, only the 
results for m(T, H) of ferromagnets [TbI ] are available. 

The 1/D expansion also works in the situations with 
nonzero magnetic field, which are not amenable to the 
methods of Refs. |l|,|l| imposing an external condition 
m = 0. An especially interesting issue is the singular be- 
havior of x{H, T) for H,T ^ for the square-lattice and 
linear-chain models. For any iJ ^ 0, the spins with low- 
ering temperature come into a position nearly perpendic- 
ular to the field, thus hm//^o limr^o xi^^: T) = 1/(2 Jq), 
where Jq is the zero Fourier transform of the exchange 
interaction, Jq = zJ , z is the number of nearest neigh- 
bors. This value coincides with the susceptibility of the 
three-dimensional classical antiferromagnets on bipartite 
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lattices in the direction transverse to the spontaneous 
magnetization. For H = 0, the spins assume all di- 
rections, including that along the infinitesimal field, for 
which the susceptibility tends to zero at T — > 0. Thus 
liinT^olimH^oxiH,T) = l/(2Jo)(l - 1/D). One can 
see that the difference between these two results is cap- 
tured exactly in the first order in 1/D. According to Rcf. 
| p2[ , for any H with lowering temperature x{H,T) 
increases, goes through the flat maximum, decreases, at- 
tains a minimum and then goes up to the limiting value 
l/(2Jo). 

The existence of the interesting features described 
above, which should be also pertinent to quantum anti- 
ferromagnets, have never been checked numerically. That 
is why we have undertaken MC simulations for classical 
AFMs in square-lattice and linear-chain geometries. 

Our systems are defined by a classical Heisenberg 
Hamiltonian 
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-H Si + Jij SiSj 



(1) 



where S is a D-component normalized vector of unit 
length (|S| = 1), H is a magnetic field and the exchange 
coupling Jij is J > for nearest neighbors and zero oth- 
erwise. The mean-field transition temperature is given 
by T^'^^ = Jo/D = zJ/D. Although there is no phase 
transition in our model, it is convenient to choose T^^^ 
as the energy scale and to introduce dimensionless tem- 
perature, magnetic field, and susceptibilities 
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h = H/Jq, Xa = JoXa, (2) 



where Xa = d{Sa) / dHa and a = x,y,z. 

In the limit D — > oo, the model Eq. (|l|) is exactly solv- 
able and equivalent to the spherical model. The solution 
includes an integral over the Brillouin zone taking into 
account spin-wave effects in a nonperturbative way. The 
latter leads to the absence of the phase transition for the 
spatial dimensionalities d < 2. 

The 1/D corrections to the spherical-model solution 
have been obtained in Refs. l9|-[T2|. They include double 
integrals over the Brillouin zone and are responsible for 
the maximum of the antiferromagnetic susceptibility at 
9^1 [|ll|. For small fields and temperatures, h,0 ^ 1, 
the field-induced magnetization m for the square-lattice 
model simplifies to 
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(3) 



which follows from Eqs. (4.9) and (2.23) of Ref. The 
log term of the above expression is responsible for the sin- 
gularity of both transverse and longitudinal (with respect 
to the field) susceptibilities. 



X± = m/h, XII = dm/dh, 



(4) 



which was mentioned above. For h = they have the 
form X = [1 - 1/D + e/D]/2, whereas for h ^ the 



limiting value at 6 = and the slope with respect to 6 
are different: x = {1 - [^/iT^D)] ln[16/(e'^/i2)]}/2. In the 
latter case, x has a minimum at 6* = 6** = tt/ ln(16//i^). 
There are corrections of order 9^ and 1/D^ to Eq. (|3). 
The latter renormalize the last, regular term in Eq. (g) 
(seeEq. (8.2) of Ref. ||ll|). The l/D^ corrections cannot, 
however, appear in the log term of Eq. (^), because this 
would violate the general properties of x{H, T) discussed 
above. 

For the linear chain, the magnetization in the region 
/i, <C 1 to first order in 1/D is given by |n2[ 
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(5) 



The transverse susceptibility of the linear chain behaves 
qualitatively similarly to that of the square lattice. The 
minimum of x± is attained at 9 = h"^^^ which is smaller 
than in two dimensions. The longitudinal susceptibil- 
ity X|| corresponding to Eq. (|^) has a minimum at 9 = 
31/3^2/3 ^ 1^ g^^^ ^ maximum at 9 = S-^/^/i^/^ -C h. 

For comparison, the zero-field Takahashi's results pl| ] 
for the Heisenberg model on the linear chain and square 
lattice can for 6* ^ 1 be rewritten in the form [O 
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where the exponentially small terms are neglected. For 
both lattices the low-temperature expansion is the same 
to order 9: x = (1/3) + (1/9)6* -I- . . ., and the results 
diverge at 6* ~ 1. The coefficient in front of 9 here is at 
variance with the l/Z?-expansion results above for D = 
3. It was argued in Ref. [ll[ that the correct general- 
form of the low-temperature expansion of the zero-field 
susceptibility for both square lattice and the linear chain 
reads 

i.e., it is reproduced to order 9 at the second order of 
the 1/D expansion. This formula is in accord with Taka- 
hashi's theory. 

In order to check the validity of the analytic results 
from the 1/D expansion above for the most realistic case 
of D = 3, we performed Monte Carlo simulation for 
three-component classical spins on a chain with length 
N as well as on a square lattice of size N ^ L x L, both 
with periodic boundary conditions. In our Monte Carlo 
procedure, a spin is chosen randomly and a trial step is 
made where the new spin direction is taken randomly 
with equal distribution on the unit sphere. This trial 
step does not depend on the initial spin direction. The 
energy change of the system is computed according Eq. 
(|l|) and is accepted with the heat-bath probability. One 
sweep through the lattice and performing the procedure 
described above once per spin (on average) is called one 
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FIG. 1. Temperature dependence of the longitudinal and 
transverse susceptibility for the square lattice for different 
values of the magnetic field h. The points are results from 
Monte Carlo simulations for L = 64 and h = 0, 0.1, and 1. 
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FIG. 2. Temperature dependence of the longitudinal and 
transverse susceptibility for the chain for different values of 
the magnetic field. The points are results from Monte Carlo 
simulations for L = 100 and ft = 0, 0.1, 0.3, and 1. 



Monte Carlo step (MCS). We start our simulation at high 
temperature and cool the system stepwise. For each tem- 
perature we wait 6000MCS (chain) and 4000MCS (square 
lattice), respectively, in order to reach equilibrium. Af- 
ter thermalization we compute thermal averages (. . .) for 
the next 8000MCS (chain) and 6000MCS (square lattice), 
respectively. 

The relevant quantities we are interested in are the 
magnetization m = jriz = {Mz) and the components of 
the susceptibihty Xa = -^((Af^) — {Ma)^), where the z 
axis is directed along H, a = x,y, z, and Ma = jfj^i^?- 
We have used the formula above for Xa to simulate the 
zero-field and longitudinal susceptibility, x\\ = Xz- For 
the transverse susceptibility, x± = Xx = Xy, nonzero 
field it is more convenient to use Eq. (|4|). For h — Q the 
transverse and longitudinal susceptibilities are identical 
and calculated as x± = X|| = iXx + Xy + Xz)/3. 

With intent to minimize the statistical error and to be 
able to compute error bars we take averages over Nr = 
100 independent Monte Carlo runs. The error bars we 
show are the mean errors of the averages a / y/N^, where 
a is the standard deviation of the distribution of thermal 
averages following from the independent runs. 



We start the comparison of theoretical and numerical 
results with the square lattice. Fig. |^ shows the tem- 
perature dependence of the reduced longitudinal suscep- 
tibility x\\ ^'^d reduced transverse susceptibility x_l for 
different values of the magnetic field, both for the sys- 
tem size L = 64. The corresponding results for the spin 
chain with system size L — 100 are presented in Fig. |2[ 

We investigated possible finite-size effects by varying 
the lattice size. However, we did not find any significant 
change of our data for lattice sizes in the range L = 
16 ... 64 (square lattice) and L = 40 . . . 100 (linear chain). 
Also, we did not find any systematic change of our results 
for longer Monte Carlo runs so that we believe to present 
data corresponding to thermal equilibrium. 

Note, that for all Monte Carlo data shown the er- 
ror bars of the transverse susceptibility are smaller than 
those of the longitudinal one since the transverse sus- 
ceptibility follows directly from the z component of the 
magnetization while the longitudinal susceptibility is cal- 
culated from the fluctuations of the z component of the 
magnetization. In the case h = the transverse and 
longitudinal susceptibility are identical and follow from 
fluctuations of the magnetization so that the error bars 
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are larger. 

For the square lattice as well as for the chain the nu- 
merical data confirm the non-analytic behavior of x in 
the limit of temperature T ^ 0, i. e. the limiting values 
X-L = X|| = 1/2 for /i 7^ and X-L = X|| = 1/3 foi' h = Q. 

Especially for the square lattice, the Monte Carlo data 
agree reasonable with the first-order 1/D expansion in 
the whole range of temperatures. On the other hand, at 
low temperatures the agreement with Takahashi's theory 
within error bars is achieved. Our numerical data thus 
confirm that the coefficient in the linear-0 term in x in 
Takahashi's theory is accurate. For h — 1 and 6* > 1, the 
MC data fall slightly below the l/D-expansion curve. 
Both are again in accord with each other for 6' > 3 (not 
shown) . 

The maximum of the longitudinal susceptibility of the 
square-lattice model for h — 1 looks much sharper than 
that of the theoretical curve. This feature, as well 
as the hump on the h = 0.1 curve at slightly lower 
temperature, are possible indications of the Berezinsky- 
Kosterlitz-Thouless (BKT) transition. The reason for 
that is an effective reduction of the number of spin com- 
ponents by one at sufficiently low temperatures in the 
magnetic field (the effect mentioned in the introduction), 
so that the Heisenberg model becomes effectively D = 2 
and it can undergo a BKT transition in two dimensions. 
We have not, however, studied this point in detail in this 
work. 

For the antiferromagnetic chain our MC simulation 
data are in a qualitative agreement with the l/D ex- 
pansion, although the discrepancies are stronger. 

Unfortunately, we could also not perform simulations 
for even lower values of the field h for the following rea- 
son: The singular behavior of x stems from the fact that 
for h > the spins tend to come into a position perpen- 
dicular to the field. For fields as small a,s h — 0.01 (curve 
4 in Figures 1 and 2) the amount of energy related to this 
ordering field is 100 times smaller than the exchange in- 
teraction energy. Therefore the corresponding relaxation 
for this energetically favorable state takes very long in a 
Monte Carlo simulation, especially for these low temper- 
atures, where this effect occurs for low fields. 

Our MC simulations showed for the first time the sin- 
gular behavior of the susceptibility of classical antiferro- 
magnets at low temperatured and magnetic fields. The 
results are in accord with predictions based on the first- 
order 1/D expansion ||ll|,|l^. It would be interesting to 
try deriving the corresponding low-temperature results 
[cf. Eqs. (|) and (|)] without using the 1/D expansion. 
One of the formulas of this type already exists: It is Eq. 
(^. A candidate among theoretical approaches is the 
chiral perturbation theory of Ref . ||^ , which is applicable 
to quantum models, as well. 

The features manifested here by classical antiferromag- 
nets should be pertinent to quantum models, as well. The 
effects observed here could be checked with the help of 
the QMC simulations which achieved recently a substan- 



tial accuracy (see, e.g., Refs. and [||). Another possi- 
bility is to map the quantum model on the classical one 
Q and to perform classical MC simulations. One should 
also mention an alternative way of mapping of quantum 
magnetic Hamiltonians on classical ones with the help of 
the coherent-state cumulant expansion ]T^ , which is a 
rigorous expansion in powers of 1/S. 
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